Isothermal-isobaric molecular dynamics using stochastic velocity rescaling 
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The authors present a new molecular dynamics algorithm for sampling the isothermal-isobaric en- 
semble. In this approach the velocities of all particles and volume degrees of freedom are rescaled by 
a properly chosen random factor. The technical aspects concerning the derivation of the integration 
scheme and the conservation laws are discussed in detail. The efficiency of the barostat is exam- 
ined in Lennard- Jones solid and liquid near the triple point and compared with the deterministic 
Nose-Hoover and the stochastic Langevin methods. In particular, the dependence of the sampling 
efficiency on the choice of the thermostat and barostat relaxation times is systematically analyzed. 



I. INTRODUCTION 

Starting from the breakthrough article of Andersen^ 
many schemes have been proposed to modify the New- 
ton equations of motion so as to perform molecular dy- 
namics (MD) in the isothermal-isobaric ensemble, where 
the number of particles, the external pressure and the 
external temperature (NPT) are fixed. Most of these 
schemes are based on an extended-Lagrangian formula- 
tion. Before proceeding to overview the details of the 
extended system method, it is better to separate two 
issues: pressure and temperature control. Within An- 
dersen's framework the control of pressure was achieved 
by adding to the Hamiltonian auxiliary dynamical vari- 
ables, which represented the volume V and its conjugate 
momentum. The resulting configurations were sampled 
from the NPH ensemble, where the enthalpy H was con- 
stant. This approach has been further extended to allow 
anisotropic variation of the cell shape and size££- Several 
almost-equivalent schemes have been proposed, with sim- 
ilar performance and similar parameters, often in combi- 
nation with different thermostats (see e.g. Refs. 0H). 

In the Andersen work, the temperature was regulated 
by additional stochastic collisions, leading to a sampling 
of the NPT ensemble. The idea underlying the Ander- 
sen's thermostat was thus similar to the use of Monte 
Carlo, and inherited from Monte Carlo the characteristic 
weak points such as discontinuity of the trajectories, lack 
of a conserved quantity, at least in the original formula- 
tion, and artifacts in the computation of dynamical prop- 
erties. Another intrinsic disadvantage of the algorithm 
due to its stochastic nature was the not-rcproducibility 
of the trajectory, unless the same random number gener- 
ator was employed. 

In subsequent refinements, the Nose-Hoover determin- 
istic thermostats^ replaced the stochastic one. For this 
purpose, a new auxiliary variable was introduced, playing 
the role of a friction on the particles and/or on the cell 
dynamics and controlling the temperature. This method 
allowed for the definition of a conserved quantity, cru- 
cial to check the integration accuracy. However, a weak 
point was that the time-evolution of the kinetic energy 
was described by a second-order differential equation, 



resulting in an algorithm which is inefficient for equili- 
bration. Even worst, the lack of ergodicity provided a 
poor control on the temperature in difficult cases such 
as harmonic solids and, in general, systems with normal 
modest Thermostat chains suggested in Ref. H partially 
solved these problems, however at the price of a more 
complex formulation and a larger number of parameters. 

A possible alternative to the deterministic thermostat 
is stochastic Langevin dynamics^ Its combination with 
extended-system control of the pressure, the so-called 
Langevin piston scheme^ was proven to be efficient and 
ergodic even in difficult cases. The lack of a conserved 
quantity has been a traditional drawback of stochastic 
MD, which was however recently solved^ However, it 
is worth noting that the choice of the friction drastically 
influences the dynamics and, indirectly, the sampling ef- 
ficiency. On the one hand, a too small friction coefficient 
may cause a very long equilibration time. On the other 
hand, a too large friction coefficient may induce strong 
perturbations of the system dynamics hindering the ex- 
ploration of the phase space. 

A further possibility for temperature and pressure con- 
trol is the deterministic scheme suggested by Bcrcndscn 
et al*£ This scheme is not based on an extended La- 
grangian, and the volume of the simulation box is driven 
by the difference between the internal pressure and the 
target one. Similarly, the kinetic energy is driven by its 
deviation from the target one. The algorithm is intu- 
itive, easy to implement and efficient in the equilibration 
phase. Unfortunately it cannot be used when the correct 
distribution is required (e.g. when fluctuations have to be 
calculated), since it does not sample the NPT ensemble. 
Moreover, in difficult cases the Berendsen scheme system- 
atically transfers energy to the slowest degrees of freedom 
and can lead to wrong results . 11 ' 14 Hence the "common 
wisdom" suggests that the initial equilibration should be 
carried out with the Berendsen scheme, while the ensem- 
ble averages should be calculated from the subsequent 
runs performed with an extended-system method. 

In a recent paper— we proposed a thermostat which 
can be regarded stochastic version of the Berend- 
sen scheme. Similarly to the Nose-Hoover scheme,S&>l 
it is a global thermostat, in the sense that it acts only 
on the total kinetic energy of the system, and produces 
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configurations in the canonical ensemble. In addition, 
this scheme allows defining a conserved quantity, does 
not suffer of ergodicity problems in solids^ and has been 
used successfully in practical applications for equilibra- 
tion purposes^ or to perform ensemble averages! 17 ! 18 Fi- 
nally, it has a minimal impact on the system dynamics, 
thus being an optimal choice for the computation of dy- 
namical properties^ 

The object of the present article is to combine the 
global stochastic thermostat introduced in Ref. [TBI with 
a standard barosta^i^ to sample the isothermal-isobaric 
ensemble. The properties of the resulting algorithm are 
discussed in details and illustrated on Lennard-Jones 
(LJ) liquid and solid systems. A special attention is 
dedicated to the parameterization of the thermostat and 
barostat relaxation times. In particular, we provide a 
framework for the choice of the barostat relaxation time 
which is based on the conservation of the total effective 
enthalpy, thus on the accuracy of the resulting sampling. 
For the thermostat relaxation time, we discuss its impact 
on the sampling efficiency by measuring the autocorre- 
lation times of several relevant physical quantities. In 
this sense, we compare the new scheme with the Nose- 
Hoover^ii and with the Langevin piston barostats^i 

This paper is organized as follows: In Section II we 
formulate the most important aspects of the theoretical 
framework linking statistical theory to equations of mo- 
tion. In particular, we demonstrate that the stochastic 
rescaling algorithm samples the isothcrmic-isobaric dis- 
tribution. Section III contains several applications of our 
scheme to LJ fluid and solid systems. Some comments on 
the practical implementation of the barostat and a reli- 
able discretization scheme are collected in the Appendix. 



II. THEORY 



center of mass of the system is conserved. Thus, it is 
convenient to adopt as dynamical variables the center- 
of-mass coordinate q cm and momentum p cm , and the co- 
ordinates r and momenta tt relative to the center of mass 



Pcm = ^ ^ Pj 



Ti = q t - q cm 

TTi = Pi 



(2a) 

(2b) 

(2c) 
(2d) 



After the change of variables, the NPT ensemble reads 



(3) 



where K(tt) = ^\ |7ri| 2 /(2m.;) is the kinetic en- 
ergy relative to the center of mass and K cm {p cm ) = 
|Pcm| 2 /(2 J2i m i) i s the kinetic energy of the center of 
mass. The delta-functions in Eq. ([3]) are related to the 
conservation of total momenta and inertia in Eqs. © 
and to the fact that, in the center-of-mass framework, 
J2j 171 j r j — and J2j = 0- As a further simplifica- 
tion, the redundant variables q cm and p cm and can be 
integrated out, resulting in the distribution 



V N pt{^^V) oc Ve -mW+u{r,v)+i 



x 6 



(4) 



A. Isothermal-isobaric ensemble 

We consider a system of N particles, described by co- 
ordinates <7i, momenta Pi, masses mj, contained in a box 
of variable volume V subject to an external pressure P ex t- 
With q and p we indicate the set of coordinates qi and pi, 
respectively. In the isotropic NPT ensemble the proba- 
bility distribution is usually defined a s 20 i 21 



V N PT(p,q,V) oc e 



-f3[K(p) + U(q,V)+P^ t V] 



(1) 



where K(p) = J2i \Pi\ 2 /{^ m i) is the kinetic energy, 
U(q,V) is the potential energy and (3 = (fc^T) -1 , with 
ks the Boltzmann constant and T the temperature. The 
expression for Vn pt has been addressed in the literature 
and is not unique. For instance, Attard 2 ^ proposed to use 
the logarithm of the volume as the integration variable, 
or, equivalently, an extra V" 1 term in the probability dis- 
tribution. Nevertheless we will adopt the more standard 
definition in Eq. |T]). 

When MD simulations are performed in absence of ex- 
ternal forces (i.e. . dU / dqj = 0) the momentum of the 



Here the additional V term comes from the fact that the 
size of the domain on which q cm is integrated is propor- 
tional to V. The contribution of the extra V factor is 
almost negligible for systems with more than a few tens 
of particles, though we will include it to have a formally 
correct ensemble. 

The distribution in Eq. ^ provides ensemble averages 
which are equivalent to those from Eq. (TT]), and should 
be used as a target ensemble for any MD simulation in 
the NPT ensemble when external forces are not present. 
In the next two Subsections we provide an algorithm to 
sample this distribution. 



B. Pressure control 

We first consider the problem of fixing thepressure, 
adopting a solution similar to that of Refs. OHE Since 
we work in the reference frame of the center of mass, 
we choose an initial configuration satisfying ^ tzi = 
and miTi = 0. Then wc evolve it according to the 
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following equations of motion 



Vi =- 



in, 

TTi =f i - Wi 

3[V{P mt -P ext ) + 2p- 1 } 



'/ 

V =Wr). 



W 



(5a) 
(5b) 
(5c) 
(5d) 



The long time distribution for the ensemble can be found 
solving the equation UP = 0. As it can be verified by 
direct substitution, the following class of functions solves 
LVnph = 

Pnph(k, r, V, rf) oc V~ 1 5 irnr^j 5 n^j 



Here /j = —dU/dri are the forces, rj is proportional to 
the relative change rate of the volume and W is the iner- 
tia of the piston, which determines the typical time scale 
for volume changes. The instantaneous internal pressure 
Pint is given by 



Pi t 



2K 



dU_ 

dV' 



(G) 



where the derivative d/dV is performed at fixed scaled 
coordinates, as discussed in more detail in the Appendix. 
Equations ([5]) are similar to those introduced by Hoover^ 
but for an additional term in Eq. (|5c[) . The original 
Hoover formulation, as observed by the author himself, 
gives a slightly wrong ensemble. With our correction, 
when the dynamics is combined with a thermostat as 
shown in the next Subsection, the exact NPT distribu- 
tion is sampled. An alternative solution to this problem 
have been proposed by Martyna et aZ£. 

It can be easily shown that the following quantities are 
conserved during the dynamics: 



5> = o 

i 



= o 



(7a) 
(7b) 



H = K{n) + U{r, V) - 2/T 1 log V + P ext V + 



Wrf 2 



Hn 



(7c) 



The first two conservation laws come from the absence of 
external forces. The last quantity H is close to the en- 
thalpy of the system (K+U +P ex tV), and its initial value 
Hq depends on the initial configuration and velocities. 

When an ensemble of systems is evolved according to 
Eqs. ([5|), their distribution evolves according to a corre- 
sponding Liouville-like equation. Defining the (6N + 2)- 
dimcnsional vector x = (r, 7r, V, 77), Eqs. (fSJ) can be writ- 
ten in compact form as x = G(x), and the Liouville-like 
equation on the density as p = —d x (Gp) = —Lp, where 
L is the Liouville operator. Written explicitly in terms 
of the phase-space variables r, tt, V and 77, the Liouville 
operator reads 



+ 



^— ' m,i OTi £ -^ / 07Vi \ 

i i \ i 

3[V(P int - Pext) + 2/r 1 ] d 



_d_ 

'dr. 



w 



dv > 3r) + 3V W (8) 



x 6 (K + U - 2/3- 1 logy + P ext V 



Wrj 2 



Ho 



(9) 



The delta functions here come from the conservation laws 
in Eqs. ([7]). Assuming that the system is ergodic an MD 
simulation based on Eqs. ([5]) will produce configurations 
according to this distribution. An alternative derivation 
of Eq. © can be done introducing the phase-space com- 
pressibility as it is done for instance in Ref. l23l 

Rigorously speaking, Eq. (J9j) is not the standard NPH 
distribution. Indeed, here the constant of motion is H, 
which, as discussed above, slightly deviates from the en- 
thalpy. Moreover, there is an additional V~ Y prefactor. 
However, since we are interested in NPT sampling, we 
ignore these discrepancies and proceed in the combina- 
tion of the NPH scheme with a thermostat. 



C. Temperature control 

We now combine the presented constant-enthalpy 
scheme with a thermostat. We first extend the NPT 
ensemble of Eq. ([4]) so as to include 77 as a dynamical 
variable. Since 77 is a velocity, we associate to it a kinetic 
energy and define the extended NPT ensemble as 



VNPT(n,r,V,r]) oc VS 



,-/3[K' (Tr,r)) + U(r,V)+P e *tV] 



(10) 



where K*(ir, rf) = K(n) + ^-j- is the total kinetic energy, 
which includes also the barostat kinetic energy. When 
the rj variable is integrated out, this distribution is iden- 
tical to the target one in Eq. Thus, we will use 
this distribution as a target ensemble for the NPT al- 
gorithm. It is easy to verify that this distribution is 
stationary with respect to the Liouville-like operator in 
Eq. ([5]), i.e. LVnpt(^, r, V, 77) = 0. The only reason why 
a time average performed with Eqs. ([5]) will not give re- 
sults consistent with the NPT ensemble is that Eqs. © 
have an extra constant of motion (namely H) . To sam- 
ple the NPT ensemble one has to introduce a thermostat, 
which changes at every step the value of H in such a way 
that the NPT ensemble is stationary. Since Vnpt is the 
product of a term depending on the potential energy and 
a term depending on kinetic energy, the thermostat can 
be designed to act only on the latter one, whose distri- 
bution in the Gibbs ensemble is known a priori and is 
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V NPT (K*) = (K*y 



Here N 



f 



the number of degrees of freedom, including the volume 
(+1) and excluding the center of mass (—3). In the fol- 
lowing, we describe how this task is accomplished by the 
stochastic velocity rescalingJ^ 

In the stochastic velocity rescaling, the application of 
the thermostat consists of a rescaling of the velocities of 
the system. We here rescale not only the momenta 7r but 
also the barostat velocity 77, by a factor a: 



77 *—arj. 



(11a) 
(lib) 



This change affects the total kinetic energy K* = K + 
^r~, which includes also the barostat kinetic energy. Fol- 
lowing Ref. [l5|, we calculate at each step the rescaling 
factor by propagating the total kinetic energy according 
to the equation 



dK*(t) 



K*(t) - K* 



dt 



The propagation is done for a half timestep At/ 2 as 
described in the Appendix, so that a 2 = K*{t + 
At/2)/K*(t). Here tt is the relaxation time of the ther- 
mostat, K* = Nj/(2/3) is the average kinetic energy, and 
dW is a Wiener noised At variance with the scheme 
in Ref. [IH, here the cell degree of freedom needs to be 
counted, and the rescaling procedure affects both the ve- 
locities of the particle and the velocity of the piston 77. 

As it has been shown previously^ the stationary 
distribution associated with the stochastic dynamics in 

Eq. (fT2"j) is the canonical one, i.e. (K*)~z 1 e~" K . It 
immediately follows that, when the equation is combined 
with the constant-enthalpy dynamics, one obtains config- 
urations distributed according to the NPT ensemble in 
Eq. dTUD- 

Our algorithm consists of an alternation of constant- 
enthalpy steps [Eq. ([5])] and a thermostat step [Eq. (fT2|) ]. 
as discussed in details in the Appendix. It is worthwhile 
noting that other schemes can be implemented in the 
same fashion, simply by changing the thermostat step. 
In particular, the Langevin piston algorithmic is easily 
implemented using the thermostat step as discussed in 
Refs. Ii~2ll25l . applying the friction and noise to the parti- 
cles and the cell velocities, and removing the velocity of 
the center of mass afterwards so as to eliminate the total 
force on the center of mass and obtain the same ensem- 
ble. Similarly, a standard Nose-Hoover scheme may be 
coupled with the constant-enthalpy integrator. 



D. Control of sampling errors 

In the practical implementation of MD one integrates 
Eqs. ([51) using a finite timestep. This invariably intro- 
duces small systematic errors in the sampled distribution. 



The traditional way to control this error is to perform a 
microcanonical simulation and to verify the conservation 
of the total energy. This approach is easily extended to 
variable cell simulations in the NPH ensemble, where 
the conservation of the enthalpy is used. An equivalent 
conserved quantity can also be defined for simulations 
performed with the Nose-Hoover scheme in the NPT en- 
semble. In a recent paper we have introduced a formalism 
which allows to define a conserved quantity, called effec- 
tive enthalpy, also for stochastic dynamics. This has been 
done for our stochastic velocity rescaling^ and has been 
later extended to Langevin dynamics ^ Here we discuss 
how this concept is generalized to an effective enthalpy 
in variable cell simulations. To this aim, we first repeat 
some of the ideas introduced in the previous papers, then 
we show how the compressibility is taken into account for 
our variable cell algorithm. 

We consider the (6N + 2)-dimensional vector x = 
(p, r, 77, V) in the phase space. In principle MD equations 
would produce a continuous trajectory x(t). However, in 
the practical implementation the time is discretized so 
that only a discrete sequence of snapshots of the system 
is available, xo,x\, . . . , which is then used to evaluate en- 
semble averages. In our scheme, each point Xi is obtained 
from the previous one Xi-\ applying alternatively Eq. ([5]) 
and the thermostat. We now consider the move from Xi 
to Xi+\ as a proposal for a Monte Carlo move. Even if 
in an MD simulation all the moves are accepted, it is in- 
structive to calculate the acceptance rate for that move. 
To enforce detailed balance relative to the target distri- 
bution Vnpt{x) 1 the acceptance rate in the Metropolis 
scheme should be set equal to 



min 1, 



M(x* <-x* +1 )P NPT (x* +1 ) 



(13) 



M(x i+ i <- Xi)P N p T (xi) 

The transition matrix Mix 1 <— a;) is the probability of 
generating x' as the next point in the sequence given that 
the present point is x. The star denotes time reversal of 
velocities, i.e. if x = (p, r, 77, V) then x* = (— p, r, —77, V). 
Clearly, Vnpt(x) = Pnpt{x*)- The detailed balance en- 
forced by Eq. (TT3"]) is not the standard one but its gener- 
alization for odd and even variables discussed in Ref. [24|. 
We then associate an effective enthalpy Hi to each point 
in the sequence in such a way that the acceptance rate 

in Eq. (TBJ) is equal to min {\,e~ p ^ i + 1 ~ e[ ^'\. By sim- 
ple substitution, the change of effective enthalpy between 
subsequent points is defined as 

M(z? <- x* i+1 )V NPT (xt 



H 



i+i 



-Ht= -p- l ln- 



(14) 



M(x i+ i <- Xi)V N p T (xi) 

We now proceed to an explicit evaluation of the incre- 
ment of effective enthalpy at each MD step. Our scheme 
is a sequence of two different operations: the thermostat 
and the time evolution of Eq. ([5]). When the thermo- 
stat is applied, either stochastic rescaling or Langevin, 
its propagation is done analytically, so that detailed bal- 
ance is exactly satisfied. This gives an acceptance equal 
to 1 and a vanishing change in the effective enthalpy. 15 
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On the other hand, when Eqs. (|5|) are propagated, the 
change in effective enthalpy is not trivial and has to be 
evaluated explicitly. The change in effective enthalpy is 
decomposed in two terms as 



V N PT{Xi) 



M(x l+1 <- Xi) 
. . (15) 

The first term, substituting the ensemble probability in 
Eq. (fit)]) and the definition of enthalpy in Eq. ((Tc"|) . is 

i Vnpt(x* +1 ) 
p m — = H(x l+1 ) - H{Xij+ 



V N PT{Xi 



p-iQnVi+i-hiVi), (16) 



where Vi is the volume at step i. The second term of 
Eq. (TT5)) is calculated reminding that the rate between the 
backward transition matrix and the forward transition 
matrix is equal to the Jacobian of the transform which, 
as shown in the appendix, is V^+i/Vi. Thus 



/rMn 



M(x* 



i+li 



-/T^lnVS+i-lnVi). (17) 



M(x i+ i <- Xi) 
Combining Eqs. (TTSJ), P^)l and (TT7|) one obtains: 
H i+1 -H i =H(x i+1 )-H(x i ), 



(18) 



so that the increment in effective enthalpy is exactly 
equal to the increment in the enthalpy. 

In conclusion, the calculation of the effective enthalpy 
is very similar to the calculation of the effective energy in 
Ref.[l5|. The effective enthalpy is obtained as a cumulated 
sum of all the increments of H in the constant-enthalpy 
dynamics. Alternatively, one can obtain the effective en- 
thalpy as the total enthalpy H minus the sum of all the 
contributions to the kinetic energy given to the thermo- 
stat. In the limit of zero time step, the effective enthalpy 
is exactly conserved. In real applications, the effective 
enthalpy is expected to exhibit random fluctuations and 
a small systematic drift. In principle, if exact sampling is 
required, the change in effective enthalpy can be used to 
perform hybrid Monte Carlo simulations,— or to reweight 
properly the obtained configurations.— In practice, one 
can just monitor its systematic drift and use it to quan- 
tify the sampling errors. If the drift is too large, one 
should decrease the integration time step. 



III. EXAMPLES 

We now test the discussed algorithm on a couple of 
model systems. We consider a system of 256 particles 
interacting through LJ potential in a cubic box with 
periodic boundary conditions. Simulations are carried 
out both in the solid fee phase and in the liquid one 
at temperature /3 _1 = 0.692 and pressure P ex t = 0, 
close to the triple point. We use everywhere LJ reduced 



units. The interactions are truncated at 2.5, and the 
force is smoothed between 2.25 and 2.5. Long-range cor- 
rections which take into account the difference between 
the true LJ potential and our truncated approximation 
are added to the total energy of the system. To this 
aim, a constant radial-distribution-function is assumed 
for distances larger than 2.25,— and the corrections to 
the internal pressure are evaluated as derivatives of the 
energy correction with respect to the cell volume. The 
dynamics is computed using an in-house code. 

For convenience, the piston mass W is defined as 
Njfl^Tp, where Tp is a time describing the time-scale 
of the volume fluctuations.— For the Nose- Hoover scheme 
discussed in Section IIII| the thermostat mass is chosen 
as Q = iVy/3 _1 r|,, where the thermostat relaxation time 
tt measures the thermostat strength. 

In the following, three kinds of thermostats are con- 
sidered: (a) our scaling procedure (stochastic rescal- 
ing), which is stochastic and global; (b) Langevin piston, 
which is stochastic and local; (c) Nose-Hoover, which is 
deterministic and global. This choice allows us to investi- 
gate separately which is the effect of a global thermostat 
versus a local one and how the stochastic nature of the 
algorithm affects its performance. We investigate the role 
played by the relevant simulation parameters, namely the 
timestep At, thermostat relaxation time tt, the baro- 
stat relaxation time Tp. For the Langevin scheme, we 
choose a friction 7 = (2tt) _1 , for the reasons discussed 
in Ref. [T^. The presented averages are obtained from 10 7 
steps simulations for each possible choice of the simula- 
tion parameters. 



Timestep, barostat mass and effective-enthalpy 
drift 



We consider here the choice of the timestep At and 
of the barostat relaxation time Tp. As it will become 
clear later, the rationale behind the choice of these two 
parameters is the same, thus we discuss them at the same 
time. 

In Fig. [1] we show the dependence of the effective- 
enthalpy drift (^|r) on At, fixing the other parameters 
at Tp = 0.5 and tt = 0.2. The customary procedure in 
MD simulations is to choose a At which is as large as 
possible, with the constraint that the drift of the total 
energy is smaller than a given threshold. In our case, an 
equivalent check can be done on the effective enthalpy in- 
troduced in Subsection III Dl or. for the Nose-Hoover, on 
its conserved quantity. The exact value of the threshold 
is rather arbitrary, but still it gives a quantitative esti- 
mate of the sampling errors given by the finite timestep, 
and still can be used to compare the relative accuracy of 
two different simulations. The two global schemes (Nose- 
Hoover and stochastic rescaling) exhibit a similar drift, 
which grows with At. Also for the Langevin scheme the 
drift increases with the timestep, but faster than in the 
global schemes. For further calculations, we choose a 
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At 

FIG. 1: Drift of the effective enthalpy {^f) plotted as a func- 
tion of the timestep At for the Lennard- Jones liquid (a) and 
solid (b) systems at the triple point, for stochastic veloc- 
ity rescaling (filled circles), Nose-Hoover (open circles) and 
Langevin (diamonds) algorithms. The thermostat relaxation 
time is tt = 0.2 and the barostat relaxation time is Tp — 0.5. 
All quantities are in Lennard- Jones reduced units. 



standard value At = 0.005, even if a slightly smaller 
timestep should be used in order to achieve the same 
accuracy with Langevin dynamics. 




FIG. 2: Drift of the effective enthalpy (^|r) plotted as a func- 
tion of the barostat relaxation time rp for the Lennard- Jones 
liquid (a) and solid (b) systems at the triple point, for stochas- 
tic velocity rescaling (filled circles), Nose-Hoover (open cir- 
cles) and Langevin (diamonds) algorithms. The thermostat 
relaxation time is tt — 0.2 and the timestep is At = 0.005. 
All quantities are in Lennard- Jones reduced units. 

We now consider the choice of rp. We recall that there 
is a relationship between the mass of the particles (and of 
the barostat) and the time-step. As an example, a four- 
fold decrease of all the masses produces exactly the same 
trajectory as a two-fold increase of the time-step, thus 
giving a two times faster exploration. In the same way, 



for a fixed value of the particles mass, a lower barostat 
mass will provide a faster volume evolution at the price 
of a lower accuracy. Thus, the choice of the barostat 
mass, W = Njfi^Tp, is based on a compromise between 
a small mass and a reasonable effective-enthalpy drift, 
and is very similar to the choice of At. In Fig. [2] we 
show the dependence of the effective-enthalpy drift on 
Tp. Since the drift induced by a too small Tp is due to 
sampling errors of a single degree of freedom (the vol- 
ume), it is partially hidden by the larger drift due to the 
inexact integration of the particles degrees of freedom. 
However, it is clear that for rp > 0.5 the drift due to 
the volume is negligible. We notice that in none of the 
presented simulations we have observed an appreciable 
decoupling between the barostat and the internal degrees 
of freedom in the case of small rp. We also computed the 
autocorrelation time of a few relevant quantities, and we 
noticed that the sampling efficiency always increases as 
Tp is decreased. Thus, the final rule turns out to be a 
Tp which is as small as possible, but with a reasonable 
effective-enthalpy conservation. We thus choose rp = 0.5 
to continue our tests. 



B. Thermostat relaxation time and statistical 
efficiency 




FIG. 3: Drift of the effective enthalpy {^f} plotted as a func- 
tion of the thermostat relaxation time tt for the Lennard- 
Jones liquid (a) and solid (b) systems at the triple point, 
for stochastic velocity rescaling (filled circles), Nose-Hoover 
(open circles) and Langevin (diamonds) algorithms. The 
timestep is At = 0.005 and the barostat relaxation time is 
Tp = 0.5. All quantities are in Lennard- Jones reduced units. 

We now proceed to a systematic study of the effects 
of the other relevant simulation parameter, namely the 
thermostat relaxation time Tp. To this aim, we fix a 
timestep At = 0.005 and a barostat relaxation time rp = 
0.5, and perform simulations for a wide range of rp. As it 
is seen in Fig. [3J the effect of rp on the effective-enthalpy 
drift is rather weak, at least for the two global schemes. 
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Indeed, for these choices of the parameters, the origin of 
the drift is not in the thermostat but in the constant- 
enthalpy step. The drift in Langcvin is larger than that 
of the other schemes. For tt we cannot use the same 
criterion that we used to optimize rp. Indeed the two 
stochastic schemes (rescaling and Langevin) are stable for 
any possible choice of tt- This is due to the fact that, in 
these two cases, the thermostat equations are linear and 
can be integrated in an analytic fashion (see Appendix). 
This is not true for the NH scheme, which turns out to 
be be unstable for t t ps 0.02, at least using the present 
integration scheme. Moreover, it is not obvious that a 
smaller tt should lead to faster sampling, as it will be 
seen in the following. 

We thus proceed in an analysis of the statistical effi- 
ciency as a function of the thermostat parameter, simi- 
lar to what has been done for the NVT ensemble.— A 
quantitative measure of the time needed for performing 
ensemble averages of a given quantity X is given by its 
autocorrelation time Tx (see e.g. Ref. [28h . Since a com- 
mon goal of NPT simulation is the calculation of aver- 
ages of quantities like the enthalpy and the volume of 
the system, we consider here the cases where X = H 
and X = V. The autocorrelation time of the total en- 
thalpy also provides a rough estimate of the time needed 
for the equilibration procedure. However, it must be no- 
ticed that to achieve a fast equilibration also the autocor- 
relation time of the fluctuations of the enthalpy should 
be short. Thus, we also calculate the correlation time 
for the fluctuations of these quantities, obtained with 
X = AH 2 = (H — {H)) 2 and X = AV 2 = (V - (V)) 2 
respectively. The first one is related to the bulk modulus 
and the second to the heat capacity at constant pressure. 
Thus, these autocorrelation times provide the statistical 
efficiency in the calculation of these important physical 
observables. 

The simulation results for a number of runs are sum- 
marized in Fig. 0] and Fig. [5] for the liquid and the solid 
respectively. In all cases, the performance of stochastic 
rescaling is equal or better than that of the local Langevin 
thermostat. For large tt, the two stochastic schemes 
have the same behavior. For small tt, the performance 
of the local Langevin scheme is worse. This happens for 
both the solid and liquid systems, and indicates that tun- 
ing the performance of Langevin thermostat is very dif- 
ficult since it is efficient only for a limited friction range. 
These problems in Langevin thermostat are related to 
the fact that it is local, thus it disturbs considerably the 
trajectory, hindering diffusion (in the liquid) or thermal- 
ization of the slow modes (in the solid). This is in agree- 
ment to what has been observed in Ref. [l^ for a liquid 
in the NVT ensemble. 

On the other hand, the Nose-Hoover thermostat is 
global and, similarly to the velocity rescaling, has a very 
small impact on the trajectory. In the present calcu- 
lations, performed at the triple point, we don't expect 
the solid to be harmonic enough to introduce ergodicity 
problems. Its performances are excellent when looking 




FIG. 4: Autocorrelation times of the average enthalpy th, 
of the enthalpy fluctuations t AH 2, of the volume ry and of 
the volume fluctuations t a1/ 2 , plotted against thermostat re- 
laxation time tt, for the Lennard- Jones liquid at the triple 
point, for stochastic velocity rescaling (filled circles), Nose- 
Hoover (open circles) and Langevin (diamonds) algorithms. 
The timestep is At = 0.005 and the barostat relaxation time 
is rp — 0.5. All quantities are in Lennard- Jones reduced 
units. 




FIG. 5: Autocorrelation times of the average enthalpy th, 
of the enthalpy fluctuations t AH 2, of the volume Tv and of 
the volume fluctuations t AV 2 , plotted against thermostat re- 
laxation time tt, for the Lennard- Jones solid at the triple 
point, for stochastic velocity rescaling (filled circles), Nose- 
Hoover (open circles) and Langevin (diamonds) algorithms. 
The timestep is At = 0.005 and the barostat relaxation time 
is Tp — 0.5. All quantities are in Lennard- Jones reduced 
units. 



at the autocorrelation time of enthalpy (th) and volume 
(ry), for both the solid and the liquid. For the solid, a 
resonance effect leads to a peak in these autocorrelation 
times for tt ~ 2. However, fluctuations of the enthalpy 
{ t ah 2 ) are well thermalized only for a limited range of 
choices for tt- This is related to the fact that the ther- 
mostat is optimally coupled only to the modes with a 
given frequency, and is not able to thermalize the others. 
The fact that the average enthalpy has a short autocorre- 
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lation time is due to cancellations in the autocorrelation 
function, and is a signature of the enthalpy oscillations 
around the correct value. The damping of these oscilla- 
tions can be rather slow if tt is not chosen properly. 

The stochastic rescaling scheme combines the stability 
of the stochastic Langevin with the efficiency of the global 
Nose-Hoover. The coupling parameter tt in the stochas- 
tic rescaling scheme has the meaning of thermostat- 
strength, and a smaller tt always guarantees a better 
efficiency. On the contrary, in the Nose-Hoover scheme 
tt also determines which modes are optimally coupled 
with the thermostat, and it is more difficult to tune. 
Finally, the stochastic rescaling is comparable with the 
Nose-Hoover in computation of averages for small tt, 
surpassing it in computation of fluctuations. 



C. Dynamical properties 

It is also interesting to check if the rescaling scheme 
preserves the dynamics. To investigate this point, we 
consider the calculation of dynamical properties. 

10 _1 E i 




FIG. 6: Diffusion coefficient plotted against thermostat relax- 
ation time tt, for the Lennard- Jones liquid at the triple point, 
for stochastic velocity rescaling (filled circles), Nose-Hoover 
(open circles) and Langevin (diamonds) algorithms. The dif- 
fusion coefficient computed from the corresponding NVE run 
is 0.03. The timestep is At — 0.005 and the barostat relax- 
ation time is rp = 0.5. All quantities are in Lennard- Jones 
reduced units. 

For the LJ liquid, we calculate the self-diffusion coeffi- 
cient D from the mean square displacement. In order to 
overcome the practical problem of computing the mean 
square displacement in a variable cell, first we calculate 
it using the scaled coordinates V _1 / 3 <jfi, then we scale 
the obtained diffusion coefficient of a factor (V) 2 / 3 to 
recover the correct value. The diffusion coefficient as a 
function of the thermostat relaxation time tt is shown 
on Fig. O It is evident that both global schemes (velocity 
rescaling and Nose-Hoover) do not affect the diffusion. In 
particular, the calculated value (D = 0.03) is in agree- 
ment with that obtained in microcanonical simulations 
in similar conditions. On the other hand, when the local 
Langevin scheme is used with a too high friction the dif- 



fusion is hindered. Remarkably, the choices of tt which 
lead to a slower diffusion correspond to the choices that 
lead to a worse sampling (see the previous Subsection). 
These time can also be compared with the typical colli- 
sion time, defined as the velocity autocorrelation time t v . 
It is easy to show that t v is related to D by t v = (3mD. 
In the limit of large tt (small friction), it turns out to 
be t v = 0.04. When the collisions with the thermal bath 
have the same frequency that the natural interparticle 
collisions, the thermostat turns out to be the bottleneck 
for diffusion. 




0.0 0.2 0.4 0.6 0.8 1.0 



t 

FIG. 7: Normalized velocity-velocity autocorrelation function 
c(t) for the Lennard- Jones solid at the triple point. The three 
overlapping solid lines, which are almost indistinguishable, are 
obtained using stochastic velocity rescaling, Nose-Hoover and 
microcanonical dynamics. The dashed line is obtained using 
the Langevin thermostat. All the thermostats are used with 
a thermostat relaxation time tt = 0.1 (for Langevin 7 = 5). 
The timestep is At = 0.005 and the barostat relaxation time is 
Tp = 0.5. All quantities are in Lennard- Jones reduced units. 

In the solid, where the diffusion coefficient is vanishing, 
we opt for the calculation of the full velocity-velocity au- 
tocorrelation function. The results are shown in Fig. 
where the autocorrelation function is plotted for all the 
considered thermostats, using a reasonable thermostat 
relaxation time tt = 0.1, and for a reference microcanon- 
ical simulation at similar conditions. From the plot it is 
clear that the global schemes do not affect significantly 
the dynamical properties. Even if these simulations are 
performed in the NPT ensemble, they can be used to 
calculate correct dynamical properties. This is not true 
for the local Langevin scheme, where the autocorrelation 
function is strongly affected by the thermostat. This is 
an additional confirmation that the impact on the dy- 
namical properties is not due to the stochastic nature 
but to the locality of the Langevin thermostat. 



IV. CONCLUSIONS 

A new method has been introduced to perform MD 
simulations at constant temperature and pressure. The 
method is stochastic but does not affect the dynamical 
properties. The concept of effective energy, which pro- 
vides a constant of motion for stochastic algorithms and 
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which was already introduced in previous papers ; 12 ! 15 
has been extended in the present scheme to an effective 
enthalpy. A systematic analysis of the role played by 
the control parameter has been shown, focused on the 
choice of the thermostat and barostat relaxation times. 
Moreover, a comparison of the new scheme with standard 
Nose- Hover and Langevin barostats has been shown. The 
new scheme appear to be simpler to use and easier to 
control. Further work is required to apply the presented 
algorithm to the case of nonequilibrium MD, where the 
Gibbs' ensemble does not hold^ 



APPENDIX A: INTEGRATION SCHEME 

At variance with Ref. [lH, we follow here a time- 
reversible scheme which is consistent with the Trot- 
ter decomposition scheme i 30 ' 31 The Trotter decom- 
position is usually adopted to decompose the Liou- 
ville (or Fokkcr-Planck) operator which describes the 
dynamics of an ensemble of realizations. For in- 
stance, if the exact solution of the Liouville equa- 
tion for a finite time increment At is p(t + At) = 
e -(L 1 +L 2 +L 3 )Atp(j : ^ j-j-g a pp r0 ximated solution is p(t + 
At) « e 4 'T e~~ L2 ^ e ~£ 3 At e -L 2 ±t e -U^ p ( t y Equiva- 
lently, one can state that if the equations of motion are 
in the form x = G\{x) + G^C^O + G%{x) an approximate 
propagator can be obtained applying in the proper or- 
der the analytical solutions of the equations x = G\{x), 
x = G2(x) and x = G^(x). A necessary condition is that 
the three decoupled equations can be solved exactly. 

We here decompose the equations of motion as a sum of 
three parts: the first is the propagation of the thermostat, 
and the others are two equations of motion whose sum 
is equivalent to the NPH dynamics. The flow of the 
integration scheme is as follows: 

1. Propagate the thermostat for a time At/2. 

2. Propagate velocities for a time At/2, according to 



r-i =0 

TTi =fi 
V =0 

3[V(P mt -P ext ) + 2p- 1 } 



n 



w 



(Ala) 
(Alb) 
(Ale) 

(Aid) 



3. Propagate positions and velocities for a time At, 
according to 



r t = h ■qr l 

rrii 

V =Wr] 
1) =0 



(A2a) 

(A2b) 
(A2c) 
(A2d) 



4. Propagate again velocities for a time At/2, as in 
step 2. 



5. Propagate again the thermostat for a time At/2, 
as in step 1. 

The sum of Eqs. (|A1[) and (|A2[) is equivalent to the origi- 
nal equations of motion ([5]) . This choice for the splitting 
is arbitrary, but it is motivated by the fact that Eqs. (|A1[) 
and (|A2j) can be solved analytically also for finite time. 

We now give explicit expressions for the propagation of 
the single stages. Stages 2 and 4 are propagated solving 
Eqs. (IAT1) : 

At i[V{P mt {t)~P ext ) + 2p- 1 ] At 

v(t + T ) = m + — 



E 



Mt)-Pi(t ) f At 
Wm, 



E 



\Ut)\ 2 1 /At s 3 



Wrrii 3 



At . . , „ . , At 

iu{t + —) = Pi{t) + Mt)—. 

The instantaneous pressure here is calculated as 



(A3a) 



(A3b) 



mt[ ) ~ W dV 



1 

w 



E 



[mm 



2 



E^'( f )' r «(*) 



dU 
dV 



(A4) 



where fij(t) is the force between particles i and j and 
T*jj is their distance, evaluated by taking into account 
the periodic boundary conditions within the minimum- 
image convention . 20 ! 32 The last term dU/dV takes into 
account any explicit dependence of the potential energy 
on the volume, such as the long range corrections i 5 -^ - 
Stage 3 is propagated solving Eqs. (|A2|) : 

njt + At) = sinh ^] EM + e ^ rt(f) (A5a) 



T](t) 



V(t + At) 



p l (t + At) 



3>)(*)At 



V(t) 



-n(t)At 



Pi(t). 



(A5b) 



(A5c) 



It is also useful to evaluate the compressibility associ- 
ated to this propagators, defined as the Jacobian of the 
corresponding transforms. The transforms in Eqs. (|A3|) 
and (|A5j) are linear, with Jacobian equal to respec- 
tively 1 and e 3v ^ At = V(t + At)/V(t). Thus, when 
a volume element in the phase space is transformed by 
Eqs. (|A3[) and (|A5[) its measure is changed by a a factor 
V(t + At) /V(t). This contribution is crucial for a proper 
evaluation of detailed-balance violations as it is done in 
Subsection III Dl 

We now describe the propagation of the thermostat 
stages (1 and 5). The thermostat can be either veloc- 
ity rescaling, Langevin or Nose-Hoover. We here simply 
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state explicitly the integration schemes for the stochas- 
tic ones, without including any derivation. The complete 
derivations can be found in the original papers. 

In the case of velocity rescaling, the propagation of the 
thermostat reads i 15 ' 19 



where 



a 2 (t) = c + 



Pi (t + At/2) = a(t) Pi (t), 
•qit + At/2) = a(t)r](t) 

(l-c^SNj.^ + R^t^K* 
N*K*(t) 



(A6a) 



(A6b) 



„ s /c(l - c)K* , , , 
v '\i N*K*(t) K 1 



and 



sign[a(t)] = sign 



R(t) 



lcN* f K*(t) 
(1 - c)K* 



(A6d) 



Here c = e At /( 2T T)^ R(t) is a Gaussian number with 
unitary variance and >!?/v*-i is the sum of JVjf — 1 inde- 
pendent, squared, Gaussian numbers. Nf = 3N — 3 + 1 
does not include the center of mass but includes the vol- 
ume degree of freedom. 

In the case of Langevin dynamics, the propagation of 
the thermostat readsi 12 i 25 



Pi (t + At/2) = c Pi (t) 



r)(t + At/2) = cr)(t) 



' rrii(l — c 2 ) 




Ri(t) (A7a) 



(A7b) 



where 



-fAt/2 



Here the Ri(t) are N independent 



vectors of normalized Gaussian numbers. After the ap- 
plication of the thermostat we remove the center-of-mass 
velocity, so as to have an ensemble exactly equivalent to 
that obtained with velocity rescaling. 
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